Phase-coexistence Simulations of Fluid Mixtures by the 
Markov Chain Monte Carlo Method Using 
Single-Particle Models 

Jun Li a , Victor M. Calo b 

a Applied Mathematics and Computational Science 
Earth and Environmental Sciences and Engineering 
King Abdullah University of Science and Technology 
Thuwal, Saudi Arabia 
b Co-Director, Center for Numerical Porous Media 
Applied Mathematics and Computational Science 
Earth and Environmental Sciences and Engineering 
King Abdullah University of Science and Technology 
Thuwal, Saudi Arabia 



Abstract 

We present a single-particle Lennard- Jones (L-J) model for C0 2 and N 2 . Sim- 
plified L-J models for other small polyatomic molecules can be obtained fol- 
lowing the methodology described herein. The phase-coexistence diagrams of 
single-component systems computed using the proposed single-particle mod- 
els for C0 2 and N 2 agree well with experimental data over a wide range 
of temperatures. These diagrams are computed using the Markov Chain 
Monte Carlo (MC) method based on the Gibbs-iWT ensemble. This good 
agreement validates the proposed simplified models. That is, with prop- 
erly selected parameters, the single-particle models have similar accuracy in 
predicting gas-phase properties as more complex, state-of-the-art molecular 
models. To further test these single-particle models, three binary mixtures 
of CH 4 , C0 2 and N 2 are studied using a Gibbs-iVPT ensemble. These re- 
sults are compared against experimental data over a wide range of pressures. 
The single-particle model has similar accuracy in the gas phase as traditional 
models although its deviation in the liquid phase is greater. The simplified 
model improves the computational efficiency significantly, particularly in the 
case of high liquid density where the acceptance rate of the particle-swap trial 
move increases. The MC method based on Gibbs- NVT ensemble is a viable 
alternative to simulate phase-coexistence of fluid mixtures. We compare, at 
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constant temperature and pressure, the Gibbs-iVPT and Gibbs-iV^T en- 
sembles to analyze their performance differences and results consistency. As 
theoretically predicted, the agreement between the simulations implies that 
Gibbs-NVT can be used to validate Gibbs-iVPT predictions when experi- 
mental data is not available. 

Keywords: single-particle model, molecular simulation, Markov Chain 
Monte Carlo method, Gibbs ensemble, phase coexistence, fluid mixtures 



1. Introduction 

The properties of phase-coexistence are important for many industrial and 
engineering applications such as the mixture separation through distillation 
column p], the transportation instability due to blockage by natural gas 
hydrates [2] or sulfur deposition [3j, CO2 sequestration and enhanced 
oil recovery [5]. To obtain these data through experimental observations 
is time consuming and expensive. Thus, molecular simulations based on 
the Monte Carlo method are auxiliary tools commonly used to understand 
phase-coexistence properties. 

The Markov Chain Monte Carlo method proposed by Metropolis et al. [6] 
is successful in simulating problems at equilibrium state and here we refer 
to it as the Monte Carlo (MC) method. It uses the importance sampling 
idea to generate configurations X, which is a high-dimensional vector made 
up of many molecular positions, according to the probability distribution 
function f(X). The consecutive configurations constitute a Markov Chain. 
The MC method estimates the expected values of the quantities of interest by 
averaging over the sampled configurations. The use of Markov chain makes 
the algorithm simple and universal but also leads to high correlation of the 
consecutive samples, which significantly increase the stochastic error in the 
MC results. Recently, the relationship of the stochastic error with the sample 
size and sampling interval was analyzed [7]. 

In MC simulations, hundreds or thousands of molecules are distributed 
inside a cubical box. Periodic boundary conditions are used to analytically 
enlarge the computational domain as it studies the behavior of a bulk fluid 
far away from the interface. For problems where the quantities of interest 
(i.e., pressure, density, mole fraction of each component) depend on molecu- 
lar position but are independent of the molecular velocity, the MC method 
records and updates only molecular positions. The MC method based on the 
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Gibbs- iVV^T ensemble was proposed in [5]. It uses two simulation boxes, one 
for liquid phase and one for the gas phase. The temperature, T, the total 
number of molecules in the two boxes, N, and the total volume of the two 
boxes, V, are fixed. The algorithm allows molecules to swap from one phase 
to the other and volume exchange between the two phases, by changing one 
box's volume and correspondingly modifying the other's volume keeping the 
total volume constant. The Gibbs- iWT ensemble MC method effectively 
simulates phase-coexistence of single component systems but becomes incon- 
venient in simulating multi-component systems as the pressure is an output 
of the simulation rather than an input parameter. For the multi-component 
systems, we use the Gibbs-iV PT ensemble MC method [H] where the pres- 
sure, p, of the two boxes is freely selected and fixed during the simulation. 
The total volume is not conserved as the volume of each simulation box 
is changed independently. Many successful applications of the MC method 
based on Gibbs-NVT and Gibbs- iVPT ensembles have been reported in the 
literature [3 [JUHIS] - 

In this paper, the phase-coexistences of binary mixtures of CH4+CO2, 
CH4+N2 and CO2+N2 are simulated using a Gibbs- NPT ensemble method. 
We study the variation with pressure of the mole fraction of each component 
in the two phases. In order to improve the efficiency of the MC simulation, 
we neglect the intramolecular structure and model CO2 and N 2 by a single 
particle as in the traditional model for CH 4 , originally proposed in [Tj5]. The 
Lennard- Jones parameters for CO2 and N2 are determined by matching the 
experimental data in [201 EI] at a temperature far away from the critical 
temperature and then used in the whole temperature range of interest. The 
single-particle modeling idea is based on the fact that the reduced equations 
of state of small molecules are similar to each other. The single-particle 
model and the selected parameters for CO2 and N2 are verified first in the 
simulations of phase-coexistence of single-component systems by comparison 
with experimental data [20j EE] over a wide range of temperatures. This com- 
parison shows that the single-particle model of C0 2 with properly selected 
parameters has similar accuracy in predicting the gas-phase properties as the 
traditional three-particle model used in [22]. To further verify the predictive 
capabilities of the single-particle model we simulate binary mixtures. As in 
the single-component case, the MC results using the simplified model agrees 
well with experimental data [231425] over a wide range of pressures. We com- 
pare the accuracy of the single-particle model against a three-particle model 
used in [26J for CO2 in the case of the binary mixture of CH4+CO2. Again, 
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the accuracy of the single-particle model is similar to that of the more com- 
plex model in the gas phase. In addition, we present the comparison between 
the Gibbs- NPT and Gibbs- iVl^T ensemble MC methods in simulating fluid 
mixture at the same temperature and pressure. This comparison shows dif- 
ference in performances between the two algorithms. While comparing the 
average results of each ensemble method shows that they are consistent with 
each other under appropriately selected conditions. 



2. Basic algorithm of the Markov Chain Monte Carlo (MC) method 

For problems of equilibrium state, the partition function of the statis- 
tical mechanics provides the formula of f(X) and the probability density 
distribution of the system's configuration X is f(X)/ j Q f(X)dX. In phase- 
coexistence problems where the quantities of interest depend only on the 
molecular position, X is a high-dimensional vector containing the positions 
of all molecules. The pressure, density, and mole fraction, which depend 
explicitly on the molecular position, are expressed as the corresponding ex- 
pected values defined by the following integral: 

= L f(X)A(X)dX 

Lf(x)dx 

where A(X) is the transient value of the quantity of interest at a particular 
configuration X of the system. As the formula of f(X) is complicated, it 
is almost impossible to get an analytical expression for (A). Traditional 
quadrature schemes are not applicable due to the large number of nodes 
required to cover the high dimensional space Q where X is defined. 

It is convenient to use the Markov Chain Mote Carlo (MC) method [6] to 
generate consecutive configurations X{ according to f{X). The MC method 
uses only f(X) rather than f Q f(X)dX. The expected value (A) is estimated 

by the average value — y^" =1 A(XA over the n sampled configurations Xj. 

ft 

The average value converges to the expected value as the sample size n grows 
infinitely. The algorithm described in [27] of the Markov Chain Monte Carlo 
method can be summarized as follows: 

1. Initialization of the configuration X: set molecular positions almost 
uniformly inside the simulation boxes; 
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2. For each cycle: 

(a) Apply trial move algorithm: the current X is randomly changed 
to X' by trial moves. The probability density of the event (X — » 
X') in the trial move is denoted by a(X — > X'). To significantly 
simplify the algorithm the following symmetric condition 

a(X X') = a(X' X) 

is required; 

(b) Apply acceptance criterion: the new configuration X' is accepted 
if Rf (random fraction uniformly distributed in [0, 1]) is less than 
the acceptance probability acciX — > X') or rejected otherwise. If 
rejected, the two consecutive configurations in the Markov Chain 
are the same. The acceptance probability is equal to 
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1 a(X> -» X)f(X') 
' a(X -> X')f(X) 



This choice is based on the detailed balance condition for the equi- 
librium state, namely 

f(X)a(X X')acc(X X') = f(X')a(X' X)acc(X' X) 
min[l B\ 

and the fact that — ——r = f3. We have that 

mm[l, p L \ 

acc(X ^X') = min[l,/(X / )//(X)] 

if the symmetric condition 

a(X X') = a(X' X) 

holds; 

3. Sample the system for the quantities A(Xj) of interest after the transi- 
tional period required to reach a state of statistical equilibrium is over. 
Samples are collected every d cycles where d is the sampling interval 
on the Markov chain. Due to the rejection of trial moves, consecutive 
samples in the Markov chain are probably identical; 

4. Stop once sufficient samples are gathered for analysis. 

A detailed analysis leading to choices of n and d that minimize the compu- 
tational requirements (memory usage and computational time) was presented 
in0. 
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3. MC algorithm based on the Gibbs-iVVT and Gibbs-iVPT en- 
sembles 



3.1. Gibbs-NVT ensemble 

As mentioned above, each molecule is modeled as a single particle and 
we refer to them simply as particles. For the description with intramolecular 
structure, the algorithms and formulas described here should be modified 
accordingly (27]. A box is employed to represent gas phase and a second 
one represents the liquid phase while different components can be found in 
a single box. Two-component systems are discussed here and the notations 
a, b are used to represent different components. The extension to cases with 
three or more components is straightforward. For the Gibbs- iVVT ensemble 
Monte Carlo method [8] introduced in [27], we have: 



f (gi,$2,V u V 2 ,N 1)a ,N lfi ,N 2)a , 



N, 



2.6 



yN^+N^ N 2 , a+ N 2 , b ^ j ( + ^ (2) 

oc 

N lia ,\N ljb \N 2ia \N 2jb \ 

where Ni ja is the particle number of the component a inside the cubic box 
1, V\ is the volume occupied by box 1, S\ is a high dimensional vector that 
contains the positions of all particles inside box 1 normalized by the box 
size (Vl) 1 / 3 {note: the subscript i is the particle index and the total particle 
number inside box 1 is Nx >a + N^f,), (3 = 1/(/cbT), k& is the Boltzmann 
constant, T is the temperature, and U\ = Ux(Si,Vi) is the total potential 
energy in box 1 estimated by the summation of pair-wise potential energies 
Uij contributed by particles i and j contained in the same box. Similar 
notation applies to the other box and other component in Eq. ([2]). The total 
volume Vtotai = Vi + V 2 , total particle numbers Ni >a + N 2 ^ a and AT 1)6 + N 2 ,b of 
each component, and temperature T are fixed in the Gibbs-iV V^T ensemble. 
The size L = V 1 ^ 3 of the simulation box is very small and the total particle 
number Ni >a + N 2a + iV 1)6 + N 2t b is usually only about one thousand due to 
limitations of computational resources. Thus, periodic boundary conditions 
are used to analytically enlarge the simulation domain. So, the contribution 
to the total potential energy U by particle's periodic images is taken into 
consideration. Taking the box 1 as an example, the following general form is 
used to express its energy summation under periodic boundary condition [27] : 

E/ 1 (^ 1 ,y 1 ) = i + ( 3 ) 

i,j,n 
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where i and j take values from 1 to (iV lja + N 1)b ) and the factor 1/2 is used 
to correct for double counting of the pair-wise contributions, = L^Sij = 
(Vi) 1 / 3 (s*j — Sj) and n is a vector of three integers from (—00,00) through 
which we can represent the contribution by the infinite particle images. For 
example, if n = (0, 0, 1), |r^- + nLi\ = |f^j + (0, 0, L\)\ which is the distance 
between particle i and one image of particle j. In the simulation, the values 
of Si instead of fi are recorded and so the normalized particle coordinates Sj 
are unchanged in the trial move of volume change. The prime over the sum 
notation means that i = j should be excluded when n = (0, 0, 0), namely we 
consider the potential energy between particle % and its infinite images but 
particle i with itself does not contribute to the potential energy. If i ^ j, 
we consider the contribution by particles i and j with n = (0, 0, 0) as well 
as the contribution by particle i and the infinite images of particle j with 
n 7^ (0,0,0). Similarly, the transient pressure p\ of box 1 at a particular 
configuration S\ is [27J: 



(N^ + N^kyT llw du\ 
yx " V 1 W x 2^ \ dr J y ' 

i,j,n 

where r = rij = \ fij + nL\\. In the case of Lennard- Jones fluid: 

12 / \ 6" 



(5) 



We specify the values of e and a for each component of a and b. If particle 
i and j belong to different components, Lorentz-Berthelot's mixing rules are 
used to compute the cross parameters 

tab = {taatbb) 1 ^ 2 

and 

(&aa + &bb) 
cr ab ~ ^ ' 

As introduced in [27], a cutoff distance r c , which is smaller than half of 
the corresponding box size, is employed to simplify the sum operation by 
limiting the number of terms with r < r c that need to be calculated explic- 
itly. Here, we use r c = 0.45L, which implies that boxes with different sizes 
have different r c . The value of r c changes after the accepted trial moves of 
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volume change as L = V 1 ^ 3 . So, the contributions by any particle i and its 
infinite images are neglected as the minimal value of their distances is L and 
larger than r c . To simplify discussion, we refer to particle j and its infinitely 
many images as the particle set of j. To compute the summation with i ^ j 
including the potential energy between particle i and the particle set of j, we 
first calculate the normalized distance between particles i and j in each coor- 
dinate axis and then get the minimal normalized distance in each coordinate 
direction while taking the infinite images of particle j into consideration. For 
example, we first compute the normalized distance As x in the x direction by 
the normalized coordinates and Sj. The periodic length of the normalized 
coordinate at all axes is 1. Thus, As x — floor(As x ) is positive and belongs 
to [0, 1), where the function floor(As x ) returns the maximum integer which 
is smaller or equal to As x . Then, the minimal normalized distance As Xjm { D 
in the x direction is As x — floor(As x ) if it is smaller than 0.5 or equal to 
1 — [As x — floor(As x )] otherwise. The minimal normalized distances As Vtm i n 
and As 2jmin in the y and z directions are computed in the same way. Now, 
the minimal normalized distance between particle i and the particle set of 
j is As min = (As^ min + As* min + As 2 zMn ) 1/2 . For any particular % and j 
{i 7^ j as i = j is neglected due to truncation) combined with all possible n 
in the summation, we only need to check a single pair-wise interaction with 
the distance equal to As min contribute to the summations in Eqs. ((3|-([4]) and 
neglect other infinitely many terms due to the truncation with r c < 0.5L. 
Namely the normalized cutoff distance s c is smaller than 0.5. According to 
the above analysis, the number of terms with r < r c in the summations is 
finite and its contribution can be computed explicitly. 

The neglected contributions with r > r c to the summations are estimated 
by tail corrections. These tail corrections for the energy and pressure sum- 
mations of box 1 are: 



tail 



k\=a ki=a 



^N lM 3 



1 ( CTkjkj 

3 V r c,i 



(6) 



and 



p tail 



Ev^ 167rN 1M N hk2 3 

k\=a k 2 =a * 



r c,l 



^kik2 
r c,l 



(7) 



The total energy and pressure are estimated by the sums of the explicit 
summations with r < r c and the tail corrections for r > r c . Note that if 
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the components a and b have the same values of e and a, the tail corrections 
degenerate to 



tail 



and 



tail 
Pi 



te{N 1>a + N 1>b ) 



16n{N ha + N lt . 
3V? 



-e a 



-e a 



-(- 

3 V r c,i 

2 / a 

3 V r c,i 



(T 

r c ,i 
a 

r c,l 



(8) 



(9) 



which are consistent with the results of the single-component system [27] . 

As the configuration X of the probability distribution function of Eq. ^ 
contains three types of independent variables which are particle coordinates, 
box volumes, and particle numbers, three kinds of trial moves are necessary: 
particle displacement, volume change, and particle swap. These satisfy the 
ergodicity condition which requires that it is possible to visit any X' G Q 
from the current X in a finite number of trial moves. After getting the total 
energy U of each box and f(X), the acceptance probability of each trial 
move can be computed. The three types of trial moves are selected with 
predetermined probabilities, which can be adjusted during the translational 
period before reaching the thermal equilibrium state. The three trial move 
algorithms used here satisfy the symmetric condition 

a(X -> X') = a(X' -> X) 

and so the acceptance probabilities are determined simply by 



acc(X -> X') 



mm 



/(*') 



In the trial move of the particle displacement, we select one box denoted 
by m from boxes 1 and 2 with equal probability and then select a particle 
denoted by i among all particles inside the box m with equal probability. 
The new normalized coordinate of the particle i is computed by 



S; 



Si + {Rh - 0.5)Ax 



where Rfi is a random fraction distributed uniformly inside [0, 1] and Ax is 
the step size of the trial move of particle displacement. We denote the new 
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total energy of the box m by U' m . This trial move satisfies the symmetric 
condition a(si — >■ sj) = — > as (Rfi — 0.5) is distributed uniformly 
inside [-0.5, 0.5] and so we have: 

acc(si ->■ s£) = min{l, exp[-/3(£/^ - U m )}} (10) 

We change Sj to if i?/2 is less than acc(sj — > s*) where H/2 is another 
uniformly distributed random fraction. After every accepted translational 
trial moves, the particle i is placed back into the box m by periodic shifting 
if its new normalized position is outside the box m, namely at least one of 
its components is outside [0, 1]. 

In the trial move of volume change, a new variable 

x = HV1/V2) = in[vy(Kotai - K)] 

is introduced [27] as the total volume Vtotai = Vi + Vi is constant in the 
Gibbs- iVV^T ensemble with 

j{x)dv 1 = f(x)v ^ Vta ^~ Vl) d x = g(x)d x 

notal 

where g(X) = /(A)Vi(Vtotai - Vi)/Vtotai- Thus, we compute a new x' by 

X / = X + (i?/3-0.5)Ay 

where AV is the step size of this trial move. Although the value range of x' 
is (—00, 00), the value of V{ is always located inside the reasonable range of 
(0, Kotai) since 

exp(x') 



V{ = V t 



total 1 



[l + exp( X ')]' 

and correspondingly, V 2 ' = Vtotai — V{. As the trial move (x — > x') satisfies 
the symmetric condition, the acceptance probability of x' is computed by 
min[l,s(X')MX)]: 



acc(x -> x) 



mm 



l\ N\ a+Ni b+X /t//\ N 2 , a +N 2 , b +1 



lf 2 ) exp [-/? {U[ + U' 2 - U x - U 2 )) 
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We change V\ and V 2 to V[ and V totai \ — V{, respectively, if i?/ 4 < acc(x — > x')- 
Si and S2 are unchanged in this move. 

In the trial move of particle swap, we select one box denoted by m from 
boxes 1 and 2 with equal probability to remove a particle i of component k. 
Simultaneously, this particle i is inserted in the other box and placed at a 
random location. The component k of particle % is selected from components 
a and b with equal probability here. Generally speaking [9], the component 
k is selected from all components with predetermined probabilities, which 
can be adjusted during the transitional period before reaching the thermal 
equilibrium state. If the particle number N m ^ of component k inside box m is 
zero, the trial move is rejected immediately and the current configuration X 
is repeated in the Markov Chain. Otherwise, we select a particle denoted by 
% among those particles of component k inside box m with equal probability. 
Taking m = 1 for instance, we remove particle i of component k from box 
1, which changes Ni t k to Ni jk — 1 and U\ to U[. Correspondingly, we create 
a particle with its coordinate Si selected randomly and uniformly in box 2, 
which changes N2,k to A^,fe + 1 and U2 to U' 2 - This trial move satisfies the 
symmetric condition and so the acceptance probability is: 

acc(N 1:k ->■ N 1>k - 1) 

= mm |1, Vl (N 2k + 1) ~ ~ 



The formula for m = 2 is similar to Eq. (12). This trial move is accepted if 
Rf 5 < acc(N hk -»■ jV 1)fc - 1). 

The step sizes Ax and Al/ are adjusted during the transitional period to 
achieve the prescribed acceptance ratios (0.5 for example) of the correspond- 
ing trial moves and fixed later to constantly satisfy the symmetric condition 
of trial moves required by the sampling process. The step size of particle 
swap is fixed at one, namely swapping one particle each time. This makes 
the acceptance ratio of particle swap fixed and usually very low if the den- 
sity of the liquid phase is very high. The low acceptance ratio increases 
the correlation degree of consecutive samples and the statistical variance of 
the MC results. The single-particle models used here effectively increase the 
acceptance ratio of particle swap trial move. 

The normalized quantities, including the normalized number density p* = 
Na 3 /V, volume V* = V/a 3 , pressure p* = po~ 3 /e, temperature T* = Tk^/e, 
and energy u* — u/e, are used in simulations to reduce the numerical error. 
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The normalization parameters e and a can be freely selected. A convenient 
selection is to set them equal to the parameters of one component of the 
mixture. 

3.2. Gibbs-NPT ensemble 

In the Gibbs-NPT ensemble [9], the total particle numbers Ni %a +N2, a and 
Ni t b + N2 t b of each component are fixed while the total volume Vtotai = V1 + V2 
is modified during the trial move of volume change. The phase-coexistence 
pressure p^ x is specified in advance for the two boxes. We denote by pa x 
the specified constant used in the following formulas to distinguish it from 
the value computed by Eq. Q which is still valid. The computed aver- 
age pressure by Eq. Q should converge to the specified value pg x . The 
original derivation of the Gibbs- NPT ensemble given in [9] is based on the 
precondition of phase-coexistence that the temperature, pressure, and chem- 
ical potential of each component are the same for the two phases. As in the 
Gibbs- NVT ensemble, three kinds of trial moves are used in the Gibbs- NPT 
ensemble: particle displacement, volume change, and particle swap. The al- 
gorithms of particle displacement and particle swap are the same as in the 
Gibbs- NVT ensemble described above. 

For the trial move of volume change, we select one box denoted by m 
from boxes 1 and 2 with equal probability to change its volume while the 
volume of the other box remains unchanged. Since V m is the only variable 
in this trial move, the results of the isothermal-isobaric ensemble is used to 
obtain the following distribution function: 

f(V m ) oc V* m ' a+Nm < b exp[-/3( P&x V m + U m )] (13) 

A new variable rj = In V m is introduced [27] . Then, 

f(V m )dV m = f(V m )V m dr] = q{V m )dr] 

where q(y m ) = f(V m )V m . We compute a new if as 

V ' = V +(Rf 6 -0.5)AV. 

Although the value range of rf is (—00, 00), the value of is always located 
inside the reasonable range of (0, 00) as V' m = exp(r/'). As the trial move 
(77 — > rj') satisfies the symmetric condition, the acceptance probability of rj' 
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is computed by min[l, q(V^) / q(V m )]: 



acc(i] — > r]') 
= min < 1, 




) 



N m ,a+N m<b +1 



exp hSftfc (K -V m )-P {U' m - U m )} 



} 

(14) 



We change V m to V' m if P/V is less than acc(77 — > rf). 

The two boxes have almost the same uniform initial state. As the volume 
of each box is changed independently in the Gibbs-iVPT ensemble, both 
boxes are prone to remain in the liquid phase which usually has a lower 
energy U than the expected gas phase. In this case, it takes a very long 
computation time for the two boxes to split into different phases which is 
the final steady state. In order to avoid such sluggish transitional period, 
we suggest discarding the trial move of volume change during the initial 
period (for example, the first 10% of the predetermined transitional period) 
such that the two boxes split quickly into two different phases via the trial 
move of particle swap. After this initial separation period, three trial moves 
are selected according to their predetermined probabilities. In addition, the 
ratio of total particle numbers of the components should be selected such 
that the mole fraction (-/V ija + N 2l a)/(Ni )a + N 2<a + iVi,b + iV 2i ;,) is between 
x a and y a , which are the steady state mole fractions of component a in the 
liquid and gas phases, respectively. This requirement also applies to the 
simulation based on the Gibbs- NVT ensemble for multicomponent systems. 
For single-component systems, where the Gibbs-iVPT ensemble is invalid, the 
simulation based on the Gibbs- NVT ensemble requires the initial density p 
to lie between the densities of the gas and liquid phases at the equilibrium 
state. 

4. Parameter determination for the single-particle model 

Usually, CH 4 is modeled as a single particle with ecH 4 /^B = 147 K and 
c"ch 4 = 3.723 x 1CT 10 m [2H] while N 2 is modeled by two atoms and CO2 by 
three atoms with fixed bond lengths and bending angle. Following [19], we 
model N2 and CO2 by single particles to improve the efficiency of the MC 
simulation and the parameters of e and a are selected appropriately to match 
the existing experimental data. In [TjJ], the parameters used in the single- 
particle model were determined according to the mean field approximation 
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to match the experimental data. As will be explained in this section, we 
advocate a simpler procedure which is easily extensible for other molecules. 

The use of single-particle model implies that the equation of state for the 
normalized quantities 

p* = Na 3 /V, p* = pa 3 /e, T* = Tk B /e 

is unique for CH 4 , N 2 and C0 2 although their parameters of e and a are 
different. Using a single-particle model for small molecules like N 2 and C0 2 
is justified by the fact that their reduced quantities p r = p/p c , p r = p/p c , 
T r = T/T c roughly satisfy the same reduced equation of state where p c , 
p c , T c are the critical values of each component. For example, the reduced 
Peng-Robinson (P-R) equation of state is: 



3.2533p r T r 4.839p 2 a(T r ) 

(15) 



Pr ' 1 - 0.25307p r 1 + 0.50614p r - 0.064044p2 



a(T r ) = 1 + (0.37464 + 1. 54226a; - 0.26992a; 2 ) (l - ^/T T 

where the acentric factor u is determined by the critical values [29]. If we 
neglect the difference due to uj between different components, the reduced 
P-R equation of state is unique. 

When the molecule is modeled by a single particle, the normalized p*, p*, 
T* satisfy a unique phase diagram and the related L-J parameters e and a are 
used to convert the normalized quantities to values with appropriate physical 
units. Fig. [T] left gives the unique phase diagram of the normalized quantities 
with comparison by the MC results in [27] . Fig. [I] right shows the converted 
results of methane using ecH 4 /^B = 147 K and cr CH4 = 3.723 x 10~ 10 m 
compared by the experimental data [30] . 



4-1. Parameter selection for C(9 2 

For C0 2 , we select eco 2 an d o"co 2 appropriately such that the converted 
results agree with the experimental data J2D] . We choose the normalized nu- 
merical results at T* — 1 (much lower than the critical T* ~ 1.35 as in Fig. [I] 
left) for converting data since the MC results deviate from experimental data 
near the critical point. At T* = 1, the normalized gas and liquid densities are 
p* = 0.029482 and p* t = 0.70111, respectively, the normalized pressures are 
p* = 0.024923 and p* = 0.024896 (not exactly the same as p* due to stochas- 
tic noise). The density ratio is p\ j p* = 23.781. While, the experimental 
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Figure 1: Transformation of the normalized MC result. 



data [20] shows that the density ratio p/, e xp/Pg,exp is 1054.84/43.662 = 24.159 
at T = 248 K and 1045.97/46.644 = 22.425 at T = 250 K. We assume that 
the variation of density ratio with temperature satisfies a linear interpolation 
and then the density ratio of experimental data at T = 248.4 K is equal to 
23.781 of the MC simulation at T* = 1. This implies that we should select 
e C o 2 /k B = T/T* = 248.4 K such that T* = 1 is converted to T = 248.4 K 
with the density ratio being closely matched. We use the density of the gas 
phase to determine another parameter <Jco 2 since the stochastic noise in the 
liquid phase is much larger than that in the gas phase. The experimental mass 
density of the gas phase is 43.662 kg/m 3 at T = 248 K and 46.644 kg/m 3 at 
250 K. So, the mass density is 44.258 kg/m 3 at T = 248.4 K by interpolation 
and the corresponding number density is 6.02 x 10 23 x 44.258 x 1000/44 m -3 
= 6.055 x 10 26 m~ 3 . As p* = Na 3 /V where N/V is the number density, we 
obtain a C02 = 3.652 x 10~ 10 m which converts the MC result p* = 0.029482 
at T* = 1 to the experimental data p 9iexp = 44.258 kg/m 3 at T = 248.4 
K. Thus, the parameters eco 2 an d o"co 2 are determined. To further justify 
this selection, we compute the MC simulation pressure of the gas phase: 
p g = p* g e C0 Jal 02 = 0.024923 x 248.4 x 1.380622 x 10- 23 /(3.652 x 10~ 10 ) 3 Pa 
= 1.755 x 10 6 Pa where we used the constant k B = 1.380622 x 10~ 23 J/K. 
The experimental pressure is 1.6746 x 10 6 Pa at T = 248 K and 1.785 x 10 6 
Pa at 250 K and so it is 1.6967 x 10 6 Pa at T = 248.4 K by interpolation, 
which is close to the pressure 1.755 x 10 6 Pa calculated by the MC method 
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with e C o 2 /k B = 248.4 K and a C o 2 = 3.652 x 1(T 10 m. 

After setting the values of eco 2 an d cr co 2 ) we perform MC simulations 
based on the Gibbs-NVT ensemble at any physical temperature of interest. 
The MC results at some particular temperatures between 216.592 K of the 
triple point and 304.1282 K of the critical point are listed in Table [T] and 
compared with the experimental data [2U] and the MC results in the liter- 
ature [22], in which CO2 is modeled by three atoms with fixed bond length 
and each atom has charge. In the elementary physical model (EPM) |22j, 
the bending angle could be fixed or flexible but the results are very close 
and deviate from the experimental data when T is close to T c . Although the 
EPM2 obtained by rescaling the parameters of the EPM is proposed in [22] 
to improve the accuracy, the temperature used in the EPM2 is inconsistent 
with the experimental value. For example, the MC results by the EPM2 at 
228 K, 258 K, 298 K agree well with the experimental data at 221 K, 250 K, 
289 K, respectively. We choose the MC results by the EPM with fixed bend- 
ing angle for the comparison of accuracy with the single-particle model used 
here, since the simulation temperature for EPM can be accurately imposed. 
Table [T] contains the pressure of gas phase of our simulation and the liquid 
pressure is neglected due to the large stochastic errors it contains. 

Table 1: Comparisons of the phase-coexistence diagrams of CO2 between 
MC results and experimental data 



T(K) 


Experimental dat 


a 1201 




MC results in 


[22] 


MC resul 


s by single-particle model 


p (MPa) 


Pi (kg/m 11 ) 


Pi (kg/m") 


p (MPa 


N (kg/m J ) 


p, (kg/m J ) 


p (MPa) 


p g (kg/m J ) 


Pi (kg/m J ) 


228 


0.82703 


21.595 


1136.34 


0.76 


19.3 


1106 


0.98264 


25.4535 


1116.52 


238 


1.1961 


31.052 


1097.05 


0.95 


23.7 


1064 


1.2938 


32.8728 


1086.55 


248 


1.6746 


43.662 


1054.84 


1.49 


37.5 


1036 


1.7573 


44.3757 


1054.83 


258 


2.2806 


60.438 


1008.71 


1.98 


19.4 


996.9 


2.2322 


56.0803 


1019.97 


268 


3.0334 


82.965 


957.04 


2.62 


66.8 


957.3 


2.8277 


71.3246 


983.056 


278 


3.9542 


114.07 


897.02 


3.44 


89.6 


909.6 


3.5732 


91.4352 


941.542 


288 


5.0688 


159.87 


822.50 


4.50 


123.2 


850.6 


4.3436 


112.850 


896.738 


298 


6.4121 


240.90 


712.77 


5.60 


164.0 


776.0 


5.4119 


149.593 


848.864 



We use the results in Table [T] to compute the relative errors for comparison 
of accuracy. The relative pressure error is defined as (p sim — p eX p)/Pexp where 
p S i m and pexp are the values of MC simulation and experiment, respectively. A 
similar definition is used for the relative density error. The comparison of the 
absolute values of the relative errors between the MC results in the literature 
(three-particle model) and current MC results (single-particle model) are 
given in Fig. [2j 
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Figure 2: Comparisons of the phase-coexistence diagrams of CO2 and the 
errors of different molecular models. 



As shown in Fig. |2j the absolute values of relative errors of the pressure 
and gas density by the single-particle model used here are smaller than those 
of the three-particle model in the temperature range from 238 K to 278 
K. This is the range where the three-particle model agrees well with the 
experimental data. The absolute value of relative error of the liquid density 
by the single-particle model is smaller than that of the three-particle model in 
the temperature range from 228 K to 258 K. Both models deviate significantly 
from experimental data when T is close to T c 304 K. 

4-2. Parameter selection for N 2 

For N 2 , we select e^ 2 and a^ 2 based on the experimental data [21]. We 
use the normalized numerical results at T* = 1 again and so p* = 0.029482, 
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p\ = 0.70111, p*Jp g = 23.781, p* = 0.024923. The density ratio P/, exp /p 9 , cxp 
of the experimental data is 25.276 at T = 98 K and 23.342 at T = 99 K. 
Thus, the experimental density ratio determined using linear interpolation at 
T = 98.77 K is equal to 23.781 of the MC simulation at T* = 1. This implies 
that we should select e^ 2 /k B = T/T* = 98.77 K such that T* = 1 is converted 
to T = 98.77 K where the density ratio is matched. The experimental density 
is given in the unit of mol/dm 3 and the gas density is 1.0466 mol/dm 3 = 
29.3048 kg/m 3 at T = 98.77 K by interpolation and so, the corresponding 
number density is 6.02 x 10 23 x 1.0466 x 1000 m" 3 = 6.3 x 10 26 m~ 3 . We 
select a Na = 3.604 x 10" 10 m which converts the MC result p* g = 0.029482 
at T* = 1 to the experimental data p 9j0xp = 29.3048 kg/m 3 at T = 98.77 K. 
According to this selection, the gas pressure in MC simulation is 

p g = p*e Na /o& a = 0.024923 x 98.77 x 1.380622 x 10~ 23 /(3.604 x 10" 10 ) 3 Pa, 

thus, p g = 0.726 x 10 6 Pa. The experimental pressure is 0.67565 x 10 6 Pa at 
T = 98 K and 0.72566 x 10 6 Pa at 99 K. Therefore, it is 0.71416 x 10 6 Pa at 
98.77 K by interpolation, which is very close to the pressure 0.726 x 10 6 Pa 
of the MC simulation with e N Jk B = 98.77 K and a Na = 3.604 x 10~ 10 m. 

We use the values of e^ 2 and <7n 2 in the MC simulations at different tem- 
peratures between 63.1526 K at the triple point and 126.19 K at the critical 
point of N2. The comparisons of our MC results by the single-particle model 
with the experimental data [21] are listed in Table [2] which contains only the 
pressure of gas phase of our simulations. The corresponding absolute values 
of the relative errors are plotted in Fig. |3j from which we can see that the 
agreement of the MC results by the single-particle model with experimental 
data is better for N 2 than CO2. 

Another interesting verification of the parameters selected here is to com- 
pare the ratios of the temperatures at the triple and critical points respec- 
tively with the ratio of e which is used for the normalization of temperature. 
The critical temperature of methane is T Cj ch 4 — 190.9 K and the temperature 
at its triple point is T t) cH 4 = 90.68 K. We have 

T t ,N 2 /T t , C H 4 = 63.1526/90.68 = 0.696, 

T c ,n 2 /T c ,ch 4 = 126.19/190.9 = 0.661 
and both of them agree well with the ratio of 

e N2 /e CH4 = 98.77/147 = 0.672. 
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Table 2: Comparisons of the phase-coexistence diagrams of N2 between MC 
results and experimental data 



T(K) 




Experimental dat 


a [21] 


MC results 


by single- 


particle model 


p (MPa) 


p g (mol/dm 3 ) 


pi (mol/dm' 5 ) 


p (MPa) p g 


(mol/dnr*) pi (mol/dm' 5 ) 


65 


0.01740 


0.03259 


30.685 


0.018585 


0.03472 


30.537 


70 


0.03854 


0.06768 


29.933 


0.042382 


0.07434 


29.772 


75 


0.07604 


0.12638 


29.153 


0.081565 


0.1351 


29.003 


80 


0.13687 


0.21737 


28.341 


0.15119 


0.2396 


28.218 


85 


0.22886 


0.35069 


27.492 


0.23808 


0.3627 


27.384 


90 


0.36046 


0.53828 


26.595 


0.37179 


0.5500 


26.497 


95 


0.54052 


0.79504 


25.640 


0.54693 


0.7923 


25.611 


100 


0.77827 


1.1409 


24.608 


0.78229 


1.1238 


24.648 


105 


1.08331 


1.6049 


23.471 


1.0660 


1.5324 


23.581 


110 


1.46581 


2.2339 


22.184 


1.4359 


2.0940 


22.425 


115 


1.93704 


3.1162 


20.658 


1.8832 


2.8450 


21.148 


120 


2.51058 


4.4653 


18.682 


2.3947 


3.8254 


19.580 



For CO2, we have 

T c ,co 2 /T c ,ch 4 = 304.1282/190.9 = 1.593 
which agrees well with the ratio of 

ecoje&u = 248.4/147 = 1.69 

but 

Tt,co 2 /T tiC R 4 = 216.592/90.68 = 2.388 

which is due to that the isothermal line at T = T t of the reduced gas-liquid 
coexistence area of CO2 is higher than that of CH4. 

5. Simulations of phase-coexistence of binary mixtures 

We use the same notation p in the following tables to represent the pres- 
sure used in experiments and the parameter used in MC simulations. The 
computed pressure by Eq. Q in MC simulations is given in the next section 
when comparing the Gibbs-iVVT and Gibbs- NPT MC simulations. 
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Figure 3: Phase-coexistence diagram of N2 and the error by the single-particle 
model when compared with the experimental data 



5.1. Gibbs-NPT ensemble MC simulation of CH4+CO2 mixture 

First, we simulate the mixture of CH4+CO2 by the Gibbs-iVPT MC 
method using the single-particle model for CO2. The temperature is fixed at 
230 K and the variations of the mole fractions of CO2 in the two phases with 
the pressure are listed in Table [3] for comparison with experimental data [21] 
and the MC results using a three-particle model for CO2 [26]. Fig. [4] plots 
the data presented in Table |3j As we can see, the MC results using the 
single-particle model of CO2 agree well with the experimental data in the 
gas phase but have larger deviation in the liquid phase (namely in xqo 2 ) 
than using a three-particle model. This is consistent with the observation in 
Fig. [2] where the prediction by the single-particle model of CO2 is worse than 
the three-particle model EPM [22J in liquid phase. 
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Table 3: Comparisons of the mole fractions in phase-coexistence of CH4+CO2 
at 230 K 



p (atm) 


Expei 


irnental data [24 J 


MC results in [26] 


MC results 


by single-particle model 


ZCO2 


2/CO2 




2/CO2 




2/CO2 


15 


0.973 


0.601 


0.973 


0.571 


0.959 


0.684 


20 


0.950 


0.475 


0.954 


0.451 


0.917 


0.526 


32 


0.885 


0.317 


0.899 


0.313 


0.796 


0.339 


40 


0.830 


0.277 


0.838 


0.260 


0.707 


0.290 


48 


0.765 


0.249 


0.796 


0.238 


0.616 


0.247 


55 


0.686 


0.236 


0.753 


0.231 


0.535 


0.226 


65 


0.528 


0.243 


0.575 


0.214 


0.382 


0.203 
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Figure 4: Comparisons of the mole fractions in phase-coexistence of 
CH 4 +C0 2 at 230 K. 

5.2. Gibbs-NPT ensemble MC simulation of CH4+N2 mixture 

We also simulate the mixture of CH4+N2 at 160 K and N2 is modeled 
by a single-particle model. The variation of the mole fraction of N2 in the 
two phases with the pressure and the corresponding experimental data [23] 
are listed in Table |4j The same results are plotted in Fig. [5] The agreement 
of our MC results with experimental data in the liquid phase is better than 
that in the simulation of the mixture of CH4+CO2 by single-particle model 
for CO2. Nevertheless, the accuracy in the gas phase does not improve in 
spite of the fact that the MC density results of pure CH 4 and N 2 agree very 
well with experimental data (see Figs. [I] and [3]). A possible explanation for 
this model behavior is that the pressure deviation of the single-particle model 
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for a pure component system increases the computed density (namely mole 
fraction) deviation as the pressure is an input parameter in the simulation 
of mixture. 

Table 4: Comparisons of the mole fractions in phase-coexistence of CH 4 +N 2 
at 160 K 



p (MPa) 


Experimental data |23J 


MC results by 


single-particle model 




2/N 2 




vn 2 


1.9913 


0.0448 


0.1742 


0.0478 


0.1521 


2.194 


0.0684 


0.2406 


0.0742 


0.2165 


2.619 


0.1205 


0.3442 


0.1438 


0.3483 


3.038 


0.1756 


0.4184 


0.2101 


0.4340 


3.395 


0.2243 


0.4657 


0.2674 


0.4926 


3.846 


0.2820 


0.5051 


0.3414 


0.5492 
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Figure 5: Comparisons of the mole fractions in phase-coexistence of CH4+N2 
at 160 K. 



5.3. Gibbs-NPT ensemble MC simulation of CO2+N2 mixture 

Finally, we simulate the mixture of CO2+N2 at 270 K. The variations of 
the mole fraction of N 2 in the two phases with pressure and the corresponding 
experimental data are listed in Table |5} The same results are plotted 
in Fig. |6j Although both CO2 and N2 are modeled with the single-particle 
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models, the agreement of our MC results with experimental data in the liquid 
phase is good, particularly in the case of high pressure. The MC results 
deviate from the experimental data in the gas phase (namely y^ 2 ) with a 
shift of about 0.05. 

Table 5: Comparisons of the mole fractions in phase-coexistence of CO2+N2 
at 270 K 



p (atm) 


Experimental data [251 


MC results by sinj 


^le-particle model 




Vn 2 




2/N 2 


37.50 


0.0108 


0.1140 


0.0150 


0.1637 


40.68 


0.0168 


0.1598 


0.0211 


0.2090 


42.25 


0.0197 


0.1783 


0.0238 


0.2274 


45.30 


0.0263 


0.2156 


0.0303 


0.2660 


50.85 


0.0368 


0.2674 


0.0413 


0.3204 


59.70 


0.0545 


0.3280 


0.0596 


0.3826 


70.00 


0.0778 


0.3770 


0.0816 


0.4317 


82.70 


0.1080 


0.4126 


0.1127 


0.4735 


91.70 


0.1319 


0.4173 


0.1343 


0.4872 


100.71 


0.1585 


0.4188 


0.1591 


0.4976 
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Figure 6: Comparisons of the mole fractions in phase-coexistence of CO2+N2 
at 270 K. 
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6. Comparison between the Gibbs-iVVT and Gibbs-iVPT ensem- 
bles 

Theoretically, both the Gibbs-iVVT and Gibbs-iVPT MC methods are 
valid for simulating mixtures. We simulate the mixture of CH 4 +N 2 at 160 K 
as an example to show differences between two methods and how their results, 
under certain conditions, are consistent. We use eN 2 and <tn 2 to normalize 
quantities including p*, V*, p*, T*, u*. To make direct comparison between 
simulation results, we first run a Gibbs-NVT MC simulation to determine the 
system's pressure. Using this pressure value as input data, we run a Gibbs- 
iVPT simulation. Similarly to the simulations of the previous section, we 
choose the total particle number N as 1024 and the normalized total volume 
V t * otal is 6000 in our Gibbs- NVT MC simulation. V* otal is selected according 
to N such that the normalized number density at the initial uniform state 
lies between the values of the liquid and gas phases in equilibrium state. We 
set 

iVl,N 2 = iV 2 ,N 2 = H8 

and 

iVi, C H 4 = iV 2iCH4 = 394 

at the initial state making the initial mole fraction of N 2 in the two boxes 
about 0.23 which is between z Nj and y^ 2 at 2.619 MPa. The final pressure 
obtained by the Gibbs- NVT MC simulation is slightly different from 2.619 
MPa because the selections of VJ.* tal , iV ijN2 , iV 2jN2 , iVi jC H 4 , iV 2 ,cH 4 are roughly 
based on the results of the Gibbs-iVPT MC simulation at pfi x = 2.619 MPa. 

In the Gibbs- NVT MC simulation, the probability for selecting the par- 
ticle displacement trial move is 0.95, 0.0009 for volume change, and 0.0491 
for particle swap after the initial short period of 1 x 10 6 cycles before which 
the trial move of volume change is avoided. We use 2 x 10 8 cycles for the 
transitional process and adjust the trial move step sizes Ax and AV every 
5 x 10 5 cycles during the transitional process such that the acceptance ratios 
of the trial moves of particle displacement and volume change approach to 
the predetermined value of 0.5. After the transitional process, the system is 
sampled every 50 cycles and 2 24 samples are collected to calculate the average 
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values. We get the average values as 



Pg 


= 2.6793MPa 




(p; = 


0.091976), 


Pi 


= 2.6664MPa 




(p? = 


0.091535), 


P9,N 2 


= 6.1215 x 10 26 nT 


-3 


(P*g,N 2 = 


0.028656), 


Pl,N 2 


= 1.82 45 x 10 27 nT 


-3 


0*,N 2 = 


0.085406), 


P<?,CH 4 


= 1.1004 x 10 27 nT 


-3 


(Pg,CH 4 = 


0.051510), 


P/,CH 4 


= 1.0195 x 10 28 nT 


-3 


(Pl,CR 4 ~ 


= 0.47726). 



As the MC results in the gas phase contain less stochastic error, we im- 
plement the Gibbs-iVPT MC simulation at pg x = 2.6793 MPa of the gas 
pressure of the above Gibbs-NVT MC simulation such that the two sim- 
ulations are comparable. The parameters setting is almost the same as in 
the Gibbs-NVT MC simulation but we slightly modified the selection prob- 
abilities of the three trial moves to 0.95, 0.0018, 0.0482 having the selection 
probability of the volume change increased twice because the volume of each 
box is changed independently in the Gibbs- NPT MC simulation. We get 



Pa 


= 2.6796MPa 




w = 


0.091986) 


Pi 


= 2.6997MPa 






0.092678) 




= 6.1528 x 10 26 nT 


-3 


(P*g,N 2 = 


0.028802) 


P(,N 2 


= 1.8 2 37 x 10 27 nT 


-3 


(P*,N 2 = 


0.085373) 


Pg,CH 4 


= 1.09 78 x 10 27 nT 


-3 


(Pg,CH 4 = 


0.051390) 


PZ,CH 4 


= 1.0204 x 10 28 nT 


-3 


(Pi,CH 4 = 


= 0.47768) 



which are very close to the results of the above Gibbs- NVT MC simulation. 
The computed gas pressure p g = 2.6796 MPa agrees very well with the 
prescribed parameter p^ x = 2.6793 MPa. 

The evolution of the normalized p* g ^ 2 , P*,n 2 > Pg,CH 4 5 P^crv P*g-> P*-> 

V* and 

V t * are given in Fig. [7] to show the comparison between the Gibbs- NVT 
and Gibbs-iVPT MC simulations. The average values agree well with each 
other but the transient results of the Gibbs-iVPT MC simulation contain 
larger stochastic error particularly in the transient volumes of the two boxes 
as they are changed independently. But, the application of the Gibbs- NVT 
MC method in the simulation of mixture is inconvenient because the pressure 
cannot be prescribed before the simulation is performed and so the study of 
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the relationship between the mole fraction and the pressure at a fixed tem- 
perature is inconvenient. Nevertheless, the Gibbs-NVT MC simulation can 
be used for the validation of Gibbs- A^PT MC simulation when experimental 
data is not available since their simulation results should be consistent with 
each other. 

7. Conclusions 

Markov chain Monte Carlo (MC) simulations of N 2 and C0 2 are per- 
formed using a single-particle model to improve the efficiency of the MC 
simulation. The corresponding Lennard- Jones (L-J) parameters are deter- 
mined according to existing experimental data. The L-J parameters for other 
small molecules can be obtained using the procedure described herein. 

The validity of the single-particle model with the selected parameters is 
verified in the simulations of systems of pure components and fluid mixtures 
by comparison with experimental data. In the pure system of N 2 , the pressure 
and the gas and liquid densities by the single-particle model agree very well 
with the experimental data over a wide range of temperatures. For C0 2 , 
the single-particle model has comparable accuracy to the traditional three- 
particle model in predicting the gas-phase properties but has larger deviation 
in the liquid phase. In the simulations of binary mixtures of CH4+C0 2 and 
CH4+N 2 , the predictions by the single-particle model are relevant for the gas 
phase although the deviation is obvious in the liquid phase. Nevertheless, 
the prediction by the single-particle model in the liquid phase becomes better 
than that in the gas phase when simulating the mixture of C0 2 +N 2 . 

The comparison between the Gibbs-NVT and Gibbs-NPT MC simula- 
tions is made in a particular case of binary mixture to show their difference in 
performance as well as the consistency of the average results at appropriate 
conditions. Although the application of the Gibbs- NVT MC simulation is in- 
convenient for systems of mixtures, it is a useful tool for checking the validity 
of Gibbs- APT MC simulation when experimental data is not available. 
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